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Computation of Linear Comparison Equations for Stability Analysis of 

Interconnected Systems 

Soumya Kundu and Marian Anghel 


Abstract —Sum-of-squares (SOS) methods have heen shown 
to he very useful in computing polynomial Lyapunov functions 
for systems of reasonably small size. However for large scale 
systems it is necessary to use a scalable alternative using vector 
Lyapunov functions. Earlier works have shown that under 
certain conditions the stability of an interconnected system can 
be studied through suitable comparison equations. However 
finding such comparison equations can be non-trivial. In this 
work we propose an SOS based systematic procedure to directly 
compute the comparison equations for interconnected system 
with polynomial dynamics. With an example of Interacting Van 
der Pol systems, we Illustrate how this facilitates a scalable and 
parallel approach to stability analysis. 

I. INTRODUCTION 

Lyapunov functions methods have long been used in 
studying stability properties of dynamical systems [1], [2]. 
Finding a Lyapunov function for a given dynamical system, 
however, is often not an easy task. Recent advances in sum- 
of-squares (SOS) methods and semi-definite programming, 
[3]-[5], have enabled algorithmic construction of polynomial 
Lyapunov functions [6], [7]. However such sum-of-squares 
based computational methods become intractable as the 
system size grows to larger than 6-8 states [8], [9]. 

It is useful to model large-scale systems in the form 
of many interacting subsystems and study the stability of 
the full interconnected system using only the subsystem 
Lyapunov functions. There are different functional forms for 
the Lyapunov function of the interconnected system, such as 
a scalar Lyapunov function expressed as a weighted sum 
of the subsystem Lyapunov functions, or applications of 
vector Lyapunov functions and comparison principles [10]- 
[13]. Particularly the formulations using vector Lyapunov 
functions are computationally very attractive because of their 
parallel structure and scalability. Based upon the results on 
comparison equations [14]-[16], the authors in [17], [18] 
introduced the concept of vector Lyapunov functions. It was 
shown that if the subsystem Lyapunov functions and the 
interactions satisfy certain conditions, then the stability of 
the interconnected system can be studied by analyzing the 
stability of a set of linear ordinary differential equations. 
However computing these comparison equations, for a given 
interconnected system, still remained a challenge. In absence 
of suitable computational tools, analytical insights were used 
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to build those comparison equations, such as the trigonomet¬ 
ric inequalities in power systems network [19]. 

In this work we use the sum-of-squares and semi-definite 
programming methods to study the stability of an inter¬ 
connected system by computing the comparison equations. 
While this approach is applicable to any generic dynam¬ 
ical system, we choose a randomly generated network of 
modifie43 Van der Pol oscillators for illustration. Each Van 
der Pol oscillator can be represented as a two-state system 
with state dynamic equations as polynomials of degree three 
[20]. The network is then decomposed into many interacting 
subsystems. Each subsystem parameters are so chosen that 
individually each subsystem is stable, when the disturbances 
from neighbors are zero. SOS based expanding interior 
algorithm [6], [7] is used to obtain estimate of region of 
attraction as sub-level sets of polynomial Lyapunov functions 
for each such subsystem. Einally SOS optimization is used to 
compute the linear comparison equation to certify stability of 
the network under disturbances. Eollowing some brief back¬ 
ground in Sec.|II]we outline the problem statement in Sec. HID 
We present the SOS-based direct approach to computing the 
comparison equations in Sec.|IV] Sec.[V]shows an application 
of comparison equations to stability analysis of a network of 
Van der Pol systems. We conclude the article in Sec. lVIl 

H. BASIC CONCEPTS AND BACKGROUND 
A. Lyapunov Stability Methods 

Let us consider the dynamical system 

x{t)=f{x{t)), f >0, xGM", /(0)=0, (1) 

with an equilibrium at the origir0, and / ; R." —>■ R" is locally 
Lipschitz. The important notions of stability are: 

Definition 1 The equilibrium point at origin is called 

1) stable in the sense of Lyapunov (i.s.L) if 

Ve>0,35>0 s.t. ||x(0)||2<5 => ||.r(f)|| 2 <e Vf, 

2) asymptotically stable if it is stable i.s.L, and 

35>0i.f. ||x(0)||2<5 =4> lim ||.r(f )||2 = 0, 

3) exponentially stable if it is asymptotically stable, and 
3h,c,5>0s.f. |]x(0)||2<5=^|[r(f)||2<ce^^'||x(0)||2 Vf 

*We choose the Van der Pol ‘oscillator’ parameters in such a way that 
these have a stable equilibrium at origin. 

^Note that this is not a restrictive assumption, since by shifting of state 
variables, the origin can always be made an equilibrium point. 


The Lyapunov stability theorem [1], [21], also called Lya¬ 
punov’s first or direct method, presents a sufficient condition 
of stability through the construction of a certain positive 
definite function. 

Theorem 1 The equilbrium point x = 0 of the dynamical 
system in (O is stable i.s.L in ^ C K.", if there exists a 
continuously differentiable positive definite function V : > 

M (henceforth referred to as Lyapunov function) such that. 


E( 0 ) = 0 , 

( 2 a) 

E(x) >0,Vx£ ^\{0}, 

( 2 b) 

V (x) > 0,Vx £ ^. 

( 2 c) 


If V satisfies —V{x) > 0,V.xG ^\{0}, then the equilibrium 
point at origin is asymptotically stable in If). Further, the 
origin is exponentially in ^ G K" if 

3a>0, s.t. —V{x)>ocV{x)yxGS!. (3) 

Here V{x) = f{x). When there exists such a function 

V (x), the region of attraction (ROA) of the stable equilibrium 


point at origin can be (conservatively) estimated as 

{x e ^ |y(x) < y'"'“ } , (4a) 

where, := argmax {x € R" |y(x) < 7} C (4b) 

The Lyapunov function can be scaled by 7 ™'^^, so that, 

,!^:={xeR''|V(.r) < 1}, (5a) 

where, y(x) = y(x)/ 7 ”““. (5b) 


Henceforth, for simplicity, we would assume, without any 
serious loss of generality, that the ROA is estimated to be 
sub-level set of y (x) = 1 . 

B. Sum-of-Squares and Positivstellensatz Theorem 

Relatively recent studies have shown that sum-of-squares 
based optimization techniques can be utilized in finding 
Lyapunov functions by restricting the search space to sum- 
of-square polynomials [ 6 ], [7], [22], [23]. Let us denote by 
R[x] the set of all polynomials in x € R". Then, 

Definition 2 A multivariate polynomial p S R[x], x S R", is 
called a sum-of-squares (SOS) if there exists h, G R[x], i G 
{1,2,... ,r}, such that p(x) = Y!i=\ Further, we denote 

the set of all SOS polynomials in x G R” by L[x]. 

Checking if p £ R[x] is an SOS is a semi-definite problem 
which can be solved with a MATLAB® toolbox SOS- 
TOOLS [3], [4] along with a semidefinite programming 
solver such as SeDuMi [5]. SOS technique can be used to 
search for polynomial Lyapunov functions, by translating (|2]i 
to equivalent SOS conditions [3], [4], [ 6 ], [24]-[27]. An 
important result from algebraic geometry called Putinar’s 
Positivstellensatz theorem [28], [29] helps in translating the 

^We will be refening to o: > 0 in as the ‘self-decay rate’. 


SOS conditions into SOS feasibili^ problems. Then the 
Putinar’s Positivestellensatz theorenO states. 


Theorem 2 Let JF = (x £ R" |mi( x) > 0,... ,Mm(x) >0} be 
a compact set, where Uj G R[x], V7 £ {l,...,m}. Suppose 


3m£R[x], so that, 


M£{c70+Lf ^o, , Vj} ,, 

& {x £ R” |m(x) > 0} is compact. 


If p{x)>0,'ixGJf, then pG O' 0 ,O’,£L[x],V_/}. 


Often for the Vi, used in this work, the existence of m(x) 
in (| 6 ll would be guaranteed [29]. 


C. Linear Comparison Principle 

Before finishing this section, let us take a look at a nice 
result on the ordinary differential equations which helps 
form the framework of stability analysis of inter-connected 
systems via vector Lyapunov functions. Noting that all the 
elements of the vector t > 0, where A = [aij] G R'”^'", 
are non-negative if and only if a/j > 0,i f j, the authors in 
[16], [17] proposed the following result; 


Lemma 1 Let A — [aif £ have only non-negative off- 
diagonal elements, i.e. aij >0, if j. Then 

v(t) <Av(t), f > 0, v£R", v(0) = vo, (7) 

implies v(t) < w(f), Vf > 0, where 

w(t)=Aw{t), f > 0, h’£R”, w(0) = v(0) = vq. ( 8) 

This result will henceforth be referred to as the ‘linear 
comparison principle’ and the differential equation in (l 8 ]l as 
the ‘comparison equation’. 

III. PROBLEM DESCRIPTION 

Eor the rest of this work, let us make the simplifying 
assumption that the dynamical system in ([T]) is in polynomial 
fornfl denoted by / £ R[x]", and that the system in O is 
(locally) asymptotically stable. 

A. Decomposed System Model 

The dynamical system in ([1]) can be expressed in the form 
of m (> 2 ) interacting, and asymptotically stable subsystems 


V/ = 1,2,... ,m, 



Xi=fi{Xi) 

Xi G M.”', X G M.” 

(9a) 




(9b) 


giif)=0: 

, Wxi G {x € M” Xj=0yj^i} 

(9c) 

where. 

II 

,...,xlfGW‘ 

(9d) 

and 

III 

Xi n X j = 0 . 

(9e) 


i=\ 


‘^For other versions of the Positivstellensatz theorem please refer to [29]. 
^Non-polynomial dynamics can be recasted into an equivalent polynomial 
form, with introduction of additional state variables and suitable equality 
constraints [7], [24], [26], [30]. 


Here Xj represents the states that belong to the 1-th subsystem 
fi G denotes the isolated subsystem dynamics, and 

gi G R[x]"' represents the neighbor interactions. 

Let us assume that the interactions can be expressed as 

Vie {1,2,...,m}, g,(x) = J^g,j(xi,xj), (10) 

where gjj G ]R[x;,Xj]"' quantifies how subsystem affects 
the dynamics of subsystem Note that (fTOl) is not a very 
restrictive assumption, since given the choice of states x, of 
the subsystem the rest of the subsystems can always be 
chosen in a way such that (fTol l holds. We denote by 

^:={1}U{;| 3{x;,X;}, s.t. g,; (x;,X;) 0} (11a) 

and X; := x/, (11b) 

je.A 

the set of indices of the subsystems in the neighborhood of 
(including the subsystem itself) and the states that belong 
to this neighborhood, respectively. 

The next step is to characterize the stability properties of 
the isolated subsystems 

Vie {1,2,...,m}, Xi=fi{xi), x,gK"'. 

by computing a polynomial Lyapunov function V) e M [x,] for 
each 1, and the corresponding estimate of the ROA as in Q. 
An SOS based expanding interior algorithm, [6], [7], is used 
to iteratively enlarge the estimate of the ROA by finding a 
‘better’ Lyapunov function at each step of the algorithm. At 
the completion of this iterative algorithm, the stability of each 
isolated subsystem (assuming no interaction) is quantified by 
its Lyapunov function Vj G K[x,], with a final estimate of the 
domain of attraction given by 

iff := {x; G R”' |V,(x,) < 1}, VI = 1,2,... ,m. (12) 

B. Stability under Interactions 
Let us define the domain 

if “ := {x G K" I x; G Iff, Vl= l,2,...,m} , (13) 

which could be interpreted as the ROA of the ‘free’ intercon¬ 
nected system (|9l), in absence of the all the interactions. The 
disturbances coming from the neighbors can be expressed 
by the subsystem Lyapunov function level-sets. While the 
equilibrium at origin corresponds to the level sets V,(0) = 
0,V1, any disturbance (or initial condition) away from this 
equilibrium would result in positive level-sets V;(x,(0)) = 
if G (0,1] for some or all of the subsystems. 

A necessary and sufficient condition of asymptotic stabil¬ 
ity (Definition [T]i can then be translated into the condition 

Vi, V,(x,(0)) = yf Vi, ^lim^V,(x,(f)) =0, (14) 

where t > 0, are solutions of the coupled dynamics 

in ®. Even though (fTTt reduces the dimensionality of the 
problem, it still remains a generally non-trivial problem. An 
attractive, and scalable, alternative approach is to construct 
a vector Lyapunov function V : R" —?> R'" 

y(x):=[yi(xi) y 2 (x 2 ) ... vMf, (i5) 


and use a comparison equation to certify if the condition 
(fl4l i holds. Restricting our focus to the linear comparison 
principle (Lemma[T]i, the aim is to seek an A = [a^] G 


and a domain ^ C such that 

y(x)< AV(x), Vx G Si G, (16a) 

where, aij>0Vi^j, (16b) 

A = [aij] is Hurwitz, and (16c) 

S is invariant under the dynamics ([TJ. (16d) 


If there exist a ‘comparison matrix’ A = [a,j] and S C 
satisfying (Hsii, then any x(0) G S would guarantee exponen¬ 
tially convergence of y(x(f)) to the origin (Lemma [U, 

3/7,00 s.t. ||y(x(t))||2<ce-'''|l^(-*(0))|l2 , Vf>0, (17) 

which also translates into exponential convergence of the 
states themselves [10]. Note that, S C 1%^, if exists, presents 
an estimate of the ROA of the full interconnected system. 


IV. COMPUTING THE COMPARISON EQUATION 
A. Traditional Approach 

In [10], [11], [13], [19], and related works, authors laid out 
a formulation of the linear comparison equation using certain 
conditions on the Lyapunov functions and the neighbor 
interactions. It was observed that if there exists a set of 
Lyapunov functions, v,-: R"' —>■ R, Vi=l,2,...,m, satisfying 
the following conditions 


Vi G {1,2,...,m}, Bfiiufia, fja > 0 such that, 

Vxi G SiC^°, fin ||x ;||2 < v,(x,) < 77,2 ||x ,||2 (18a) 

and (Vv,y/i-<-77,3 ||x,||2 (18b) 

and if the interaction terms in (ITOl i satisfy 

Vi G {1,... , 777 }, V)G^\{/} , 3^,'2>0 such that, 

VxiGSi,VxjGSj, iVvi)^gij ^<Qj\\xj\\^, (19) 

then the following comparison equation can be formed. 


Vx{t)GS, v(x(f)) < Av(x(7)), A = [ 5 , 2 ] G R'”^'" (20a) 

where, v(x) = [vi(xi) V 2 (x 2 ) ... Vm{xm)f, (20b) 

^ = {x G R”| Xi G ViG{ 1 ,...,777}}C 1%^, (20c) 
(-Jle/fla, j = i 

and dij=< Cij/BjU jGSf\{i}, Vi,Vj (20d) 

[ 0 , 

If the ‘comparison matrix’ A = [dij] is Hurwitz, then any 
invariant domain 1% C S provides an estimate of a region of 
exponential stability of the full system [11], [19]. 


B. Motivation for Direct Approach 

While this approach provides very useful analytical in¬ 
sights into the construction of the comparison matrix A = 
\dij\, it has certain computational issues. This requires finding 
the bounds in (ITsT l and ( [T9] l, and also the Lyapunov functions 
Vi, Vi, that satisfy those. Clearly the polynomial Lyapunov 
functions, ViVi, cannot satisfy the linear bounds in (ITsT l. 






Assuming that the polynomial Lyapunov functions, V, Vi, 
we found using the expanding interior algorithm (Sec. lIII-Al) 
are quadratic, we can define v, := which would 

satisfy the conditions in ( fTSl l [11], [19]. In such a case, one 


needs to find the following bounds. 


V/,VJ G {i} ,Va; G Si,yxj G Sj , 

mWxiWi < riaWxiWl. 

(21a) 

'^ylfi < -bb WxiWl, 

(21b) 

and 2 < C6'll-*<il2 U 2 ’ 

(21c) 


for some positive scalars I7n,Ci;- Then using simple 
algebra the bounds in (ITsT l and ( fT9] l can be obtained as 


Vi, Vy e ^ {i}, fjn = y/riK, fja = 

Vb 


riB = 




and Cij = 


Cu 


'^y/W\ 


(22a) 

(22b) 


Thus the computation of each element of the comparison 
matrix A in (i20l) requires multiple optimization steps. 

We may also note that some of the bounds in and 
while convenient for analytical insights, need not be optimal 
for computing a Hurwitz comparison matrix. For example, in 
II Vvfg,; ||2 is function of both Xi and Xj but is bounded 
by using only the norm on Xj. 


C. SOS Based Direct Approach 


In this work, we are interested in Slt^i, of the form. 


V/, := 

and. Si : = 




, X € . 




(23a) 

(23b) 


Note that we exclude the boundary of the isolated subsystem 
ROA, 7^ = 1V /, for reasons explained later. The comparison 
equation in (I16al) can then be translated into 


Vi, Vi{xi) < aijVj{xj), Mx e 
i=i 

i-e., Vi< Y when Vj < V; e 


(24a) 

(24b) 


since we know that aij = 0 VJ ^ Using the Positivstel- 
lensatz theorem (Theorem |2|, with m, := {'}f — Vi{xi)) and 
= S, we can cast (i24l l into an SOS feasibility problem, 

- if, + g,) + Y {^ijVj - ^ij {fj - Vj)) e (25a) 

with (J,j G E[xi], V/G {l,2,...,.m},VjG (25b) 


where x, was defined in (fTTI i. The goal is to find the ‘optimal’ 
scalars ajj\/ij G ^satisfying ( i25l l so as to obtain the tightest 
possible bound in (I16al) . We can thus formulate the following 
SOS optimization problem. 


We propose to use SOS methods to directly compute 
the comparison equation in (IThl l. in a decentralized way by 
calculating each row of A = [aij] directly at each subsystem 
level. Note that, in (HSl), we will be using quadratic (or, 
in general, polynomial) Lyapunov functions which do not 
satisfy the bounds ([T8l) -([T9b. But we may observe that. 


Lemma 2 If there exist Lyapunov functions v,-, for each i G 
{1,2,..., m }, satisfying the comparison equation in (I20al l- 
(I20bl) . for some matrix A = [dij] with 

m 

dij >0 \/i J, and Y ^ij < 0 Vi, 
j=i 

then there exists another matrix A = [atj] satisfying the 
comparison equation (I16al) with y := V?, V i, with 

m m 

aij > 0 V/ 7 ^ J, and Y ^ij < Y < 0 V/. 
j=i j=i 

Proof: Please refer to Appendix 0 ■ 

It will be useful to note here that, an application of Gersh- 
gorin’s Circle theorem [31] says that if a matrix with negative 
diagonal elements is strictly diagonally dominant then the 
matrix is Hurwitz. We are now in a position to outline the 
SOS based procedure to directly compute the matrix A = [aij\ 
in the comparison equation (IThl l. 

= [a,j] is strictly diagonally dominant if Ljj^i < kill A'- 


V/G (1,2,...,m} , min ^ aij, subject to (i25T l . (26) 

je^i 


This simple SOS formulation helps us find the comparison 
equation (IThl l in a decentralized way, by computing each row 
of the comparison matrix A = [a,j] in a single optimization 
problem at each subsystem level. The optimization problem 
can be easily implemented on a parallel platform, with the 
complexity of the problem essentially dependent on the size 
of the largest neighborhood 

Further note that, if the minimal values of all the row-sums 
in (i26l l are negative, then the matrix A = [aij] thus found is a 
strictly diagonally dominant with negative diagonal entries, 
and hence, Hurwitz [31]. However, if Yfj=i‘^ij > 0 for any i, 
then the eigenvalues of A = [a,;] need to be computed. 

Finally a note on invariance of the domain S, which along 
with the presence of Hurwitz A = [aij] guarantees that ^ is a 
domain of exponential stability, as noted in (IThT l. According 
to [11], an estimate of the ROA can be given by 



max 



< min 



where, Pi > 0, V/G (1,2,... ,m} , 


andAp<0, p := {pi,p 2 ,. ■ ■ ,Pmf ■ 


(27a) 

(27b) 

(27c) 


Then it is easy to see that. 









V. NUMERICAL EXAMPLE 

A. Model Description 

We consider a network of nine Van der Pol ‘oscillators’ 
[20], with parameters of each oscillator chosen to make them 
individually stable. Each Van der Pol oscillator constitutes a 
subsystem, with the interconnections shown below 

{1,2,5,9} Mi: {2,1,3} Mi : {3,2,8} 

M}: {4,6,7} Mi: {5,1,6} Mi : {6,4,5} (29) 

Mi:{7,4,8,9} Mi : {8,3,7} Mi : {9,1,7}. 

The dynamics of each oscillator, in presence of the neighbor 
interactions, is given by 

V;e{l,2,...,9}, 

Xj.i=Xj 2 (30a) 

Xj,2 = -Xjj+Xjj ^ (30b) 

ke^j\{j} 

where are chosen randomly from (-3,-1) and 

the interaction coefficients , V},VkGM}\{}}, are chosen 
randomly from (—0.4,0.4). 

Using the expanding interior algorithm, we find estimates 
of the ROAs of the isolated, or ‘free’, subsystems via 
quadratic Lyapunov functions. As an example, Eig.[T] shows 
a comparison of the true ROA of the isolated subsystem 9 
and an estimate using a quadratic Lyapunov function, 

^9={{x9,UX9,2)\V9<1}, (31a) 

where, Vg = 0.595Xgj+0.227x9p.X9,2+ 0.520x9 2■ (31b) 



Fig. 1. Compaiison of estimated and true ROA for isolated subsystem 9. 


B. Exponential Stability of Isolated Subsystems 

Existence of a comparison matrix A = [a,;] requires that 
the diagonal entries an are negative, which necessitates that 

V/, 3a,' > 0, so that < — a,'V,', Vx,- G (32) 

where , V/, were defined in (l23t . Note that the condition 
(|3^ is a sufficient condition of exponential stability for the 
isolated subsystems, as in (|2l. We can use SOS optimization. 


similar to (l25T l- (l26l l. to find the maximal a,, V/, the ‘self¬ 
decay rates’, for a set of given if ,V/. Higher values of 
a,' indicates better chance of finding a Hurwitz comparison 
matrix. In Eig.|2] we show the variations of a, for each i, 
when the initial level set is varied from 0 to 1. Eor each 
subsystem, as approaches 1, a, approaches 0. This shows 
that it is not possible to obtain a Hurwitz comparison matrix 
when the initial conditions lie close to the boundary of the 
estimated ROAs, and hence the exclusion of = 1 in (|2^ . 



Fig. 2. Evolution of self-decay rates against varying initial level-sets. 


C. Comparison Equation 

We recall that two sufficient criteria for a domain & to 
be an estimate of the ROA, are that the comparison matrix 
in ( [T6l l is Hurwitz and ^ is an invariant domain under the 
dynamics (|9l). To compare the performance of the traditional 
approach and the direct approach, we need to monitor how 
well the above mentioned criteria are satisfied for a set of 
arbitrarily chosen 

While this would require an exhaustive simulation over all 
possible domains & defined in (|2^ . we choose to examine 
only those St where '^ = '^ = ‘“ = '^ = y*, i.e. 


^:= 





(33) 


for some 7*G(0, 1). Eor each 7*, and domain we compute 
the comparison matrices using the traditional and the direct 
approach. Denoting by Re (A) the real parts of the eigenval¬ 
ues of a matrix, we note that if the maximum of Re (A) is 
negative, then the matrix is Hurwitz. Eurther, by applying 
(l28l l. the domain (l3^ is guaranteed to be invariant if the 
maximum row-sum of the comparison matrix is negative. 
Fig E] shows an evolution of these two properties (maximum 
Re(A) and maximum row-sum) for the comparison matrices, 
computed using the two approaches, for a range of 7*. 

We note that both the maximum row-sum and the max¬ 
imum Re(A) generally increases as 7* increases from 0 to 
1, indicating that as the domain ^ ‘expands’, it becomes 
more difficult to certify stability. We also note that, for both 
approaches, the maximum row-sum becomes positive before 


























Fig. 3. Evolution of the properties of comparison matrices, computed using 
traditional and direct approach, against varying initial level-sets. 

maximum Re (A), indicating that the ‘invariance’ criterion is 
lost before the ‘Hurwitz’ criterion. Significantly, we also note 
that both the Hurwitz and invariance criteria are satisfied for 
a wider range of 7 * in case of the direct approach than in 
the case of the traditional one. Thus, with regards to both the 
criteria, the direct approach is seen to perform better than the 
traditional approach. 

D. Test Case 

Let us illustrate how this method can be used to certify 
exponential convergence of a given initial condition to the 
origin. The system dynamics is evolved against a randomly 
generated initial condition, and is found to be converging to 
the origin. Fig.|4] shows the evolution of the states belonging 
to subsystems - 2,6,7 and 8 . The initial condition yields the 



Fig. 4. Evolution of subsystem states under an ai'bitrary disturbance. 

following level sets, 

f = [0.08,0.56,0.58,0.31,0.08,0.61,0.18,0.45,0.14]^ 

which is then used to define the domain St, in (|2^ . Then 
the SOS-based direct approach is used to compute the com¬ 
parison matrix, A € with maximum Re(A) as -0.078, 


and Ay® < 0. The solution, w(t) £ R^, of the corresponding 
comparison equation w = Aw, w(0) = y**, is plotted against 
the actual Lyapunov level sets in Fig.|5] for subsystems 2, 
6 , 7 and 8 . The trajectories w(f) exponentially converge to 



t 


Fig. 5. Comparison of the actual Lyapunov level sets and their upper 
bounds from the linear comparison equation. 

zero and, from Lemma[T] provide an upper bound on the 
corresponding subsystem Lyapunov function level sets. 

When the same procedure is done with the traditional 
approach, we obtain a Hurwitz comparison matrix. A, with 
maximum Re(A) as -0.001, but with A(y‘')'^^ > 0, thus 
violating the invariance condition. 

VI. CONCLUSIONS AND FUTURE WORKS 

A. Conclusions 

We have presented an SOS based direct approach to 
compute the linear comparison principle for stability anal¬ 
ysis of interconnected systems. We have also discussed the 
traditional approach to obtaining the comparison equations, 
and shown how the direct approach can yield ‘better’, or 
less conservative, certificates of exponential stability. Using 
a network of Van der Pol systems we have presented a 
comparison of the two approaches. The proposed approach 
can be implemented on a suitable parallel platform where 
each row of the comparison matrix, corresponding to each 
subsystem, is computed in parallel. 

B. Future Works 

A decentralized control framework can be visualized 
where each subsystem computes a local control law that 
will guarantee satisfaction of the Hurwitz and invariance 
conditions. SOS methods can be used to extend the stability 
analysis to higher order, and more general, comparison 
equations. Also, it would be interesting to see how the use 
of higher order (for example, quartic) Lyapunov functions in 
the comparison equation affects the conservativeness of the 
stability certificates. 
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Appendix 

A. Proof of Lemma \2\ 

Since Vi = vj V; and dij > 0 \/if j, we have 

m 

Vi e {1,2,... ,m} , Vi < 2vi E aijVj 

j=i 

< IdiiVi + lYaijViVj 

j^i 

<2diiVi + Yaij(.Vi + Vj) 

jfi 

= [ d,j- + E ^ij ) + E 

V ;=i / j^i 

Choosing an = (da + T!j'=i ^ifj V i, and 0,7 = dij V/ f j, and 

recalling that Y.j=i ^ij < 0 we may conclude the proof. 


